#ifndef INT_H
#define INT_H

#include <functional>
#include <vector>
#include <queue>
#include <map>
#include <cmath>
#include <algorithm>
#include <assert.h>
#include <limits>
#include <ostream>

#include "MoreFunctional.h"
#include "SeriesLevin.h"
#include "FFT.h"
#include "IntGKPData10.h"

namespace integ {

#ifdef NDEBUG
#define integ_dbg1 if (false) (*dbgout)
#define integ_dbg2 if (false) (*reg.dbgout)
#define integ_dbg3 if (false) (*tempreg.dbgout)
#else
#define integ_dbg1 if (dbgout) (*dbgout)
#define integ_dbg2 if (reg.dbgout) (*reg.dbgout)
#define integ_dbg3 if (tempreg.dbgout) (*tempreg.dbgout)
#endif

    template <class T> 
    T norm(const T& x) { return x*x; }
    using std::norm;
    template <class T> 
    T real(const T& x) { return x; }
    using std::real;

    //#define COUNTFEVAL  // If defined, then count the number of function evaluations
#ifdef COUNTFEVAL
    int nfeval = 0;
#endif

    enum INTALG {GKP, GKPOSC, OSCINFTWOPI, OSCINFONE};
#define DEFALG GKP
    // GKP = adaptive Gauss-Kronrod-Patterson with bisection for each split
    // GKPOSC = adaptive GKP, with splits at zeros of function
    // OSCINF* = integration of oscillatory functions to infinity
    //       TWOPI = period of oscillation is 2 PI
    //       ONE   = period of oscillation is 1

    template <class T>
    struct IntRegion 
    {
    public:
        IntRegion(const T _a, const T _b, INTALG _alg=DEFALG,
                  std::ostream* _dbgout=0) :
            a(_a), b(_b), error(0.), area(0), alg(_alg), dbgout(_dbgout) {}
        bool operator<(const IntRegion<T>& r2) const {return error < r2.error;}
        bool operator>(const IntRegion<T>& r2) const {return error > r2.error;}
        void SubDivide(std::vector<IntRegion<T> >* children) {
            assert(children->size() == 0);
            if (splitpoints.size() == 0) Bisect();
            if (splitpoints.size() > 1) {
                std::sort(splitpoints.begin(),splitpoints.end());
            }
#if 0
            if (a>splitpoints[0] || b < splitpoints.back()) {
                std::cerr<<"a,b = "<<a<<','<<b<<std::endl;
                std::cerr<<"splitpoints = ";
                for(size_t i=0;i<splitpoints.size();i++) {
                    std::cerr<<splitpoints[i]<<"  "; 
                }
                std::cerr<<std::endl;
            }
#endif
            assert(splitpoints[0] >= a);
            assert(splitpoints.back() <= b);
            children->push_back(IntRegion<T>(a,splitpoints[0],alg,dbgout));
            for(uint i=1;i<splitpoints.size();i++) {
                children->push_back(
                    IntRegion<T>(splitpoints[i-1],splitpoints[i],alg,dbgout));
            }
            children->push_back(IntRegion<T>(splitpoints.back(),b,alg,dbgout));
        }
        void Bisect() {splitpoints.push_back((a+b)/2.);}
        void AddSplit(const T x) {splitpoints.push_back(x);}
        uint NSplit() const {return splitpoints.size();}

        void FindOscSplits(const std::map<T,T>& fx)
        {
            T lastzero = std::min(Left(),Right());
            T minwidth = fabs(Right()-Left())/20.;
            typedef typename std::map<T,T>::const_iterator mapit;
            mapit prevpair = fx.begin();
            mapit thispair = prevpair; thispair++;
            while (thispair != fx.end()) {
                const T& f1 = real(prevpair->second);
                const T& f2 = real(thispair->second);
                const T& x1 = prevpair->first;
                const T& x2 = thispair->first;
                if (f1 < 0. && f2 >= 0. && x2 > lastzero + minwidth) {
                    T zero = (f2*x1 - f1*x2)/(f2-f1);
                    if (zero - lastzero > minwidth) {
                        AddSplit(zero);
                        lastzero = zero;
                    }
                }
                prevpair = thispair;
                thispair++;
            }
            if (splitpoints.size() == 0) {
                integ_dbg1<<"No zeros found... bisect\n";
                Bisect();
            } else {
                integ_dbg1<<"Found "<<splitpoints.size()<<" zeros\n";
            }
        }

        const T& Left() const {return a;}
        const T& Right() const {return b;}
        const T& Err() const {return error;}
        const T& Area() const {return area;}
        INTALG Alg() const {return alg;}
        void SetArea(const T& a, const T& e) {area = a; error = e;}
        void SetAlg(INTALG _alg) {alg = _alg;}

    private:
        T a,b,error,area;
        std::vector<T> splitpoints;
        INTALG alg;

    public:
        std::ostream* dbgout;
    };

    template <class T> 
    T Epsilon() 
    { return std::numeric_limits<T>::epsilon(); }
    template <class T> 
    T MinRep() 
    { return std::numeric_limits<T>::min(); }

#if 0
    template <> 
    double Epsilon<double>()
    { return 2.2204460492503131e-16; }
    template <> 
    double MinRep<double>()
    { return 2.2250738585072014e-308; }
#endif

#ifdef EXTRA_PREC_H
    template <> 
    Quad Epsilon<Quad>()
    { return 3.08148791094e-33; }
    template <> 
    Quad MinRep<Quad>()
    { return 2.2250738585072014e-308; }
#endif

    template <class T>
    inline T rescale_error (T err, const T& resabs, const T& resasc)
    {
        if (resasc != 0. && err != 0.) {
            const T scale = (200. * err / resasc);
            if (scale < 1.) err = resasc * scale * sqrt(scale) ;
            else err = resasc ;
        }
        if (resabs > MinRep<T>() / (50. * Epsilon<T>())) {
            const T min_err = 50. * Epsilon<T>() * resabs;
            if (min_err > err) err = min_err;
        }
        return err ;
    }

    template <class UF>
    int IntGKPNA(
        const UF& func, IntRegion<typename UF::result_type>& reg,
        const typename UF::result_type epsabs, 
        const typename UF::result_type epsrel,
        std::map<typename UF::result_type,
        typename UF::result_type>* retfxmap=0)
    {
        // A non-adaptive integration of the function f over the region reg.
        // The algorithm computes first a Gaussian quadrature value
        // then successive Kronrod/Patterson extensions to this result.
        // The functions terminates when the difference between successive
        // approximations (rescaled according to rescale_error) is less than 
        // either epsabs or epsrel * I, where I is the latest estimate of the 
        // integral.
        // The order of the Gauss/Kronron/Patterson scheme is determined
        // by which file is included above.  Currently schemes starting 
        // with order 1 and order 10 are calculated.  There seems to be 
        // little practical difference in the integration times using 
        // the two schemes, so I haven't bothered to calculate any more.
        if (retfxmap) retfxmap->clear();
        typedef typename UF::result_type T;
        const T a = reg.Left();
        const T b = reg.Right();

        const T half_length =  0.5 * (b - a);
        const T abs_half_length = fabs (half_length);
        const T center = 0.5 * (b + a);
        const T f_center = func(center);
#ifdef COUNTFEVAL
        nfeval++;
#endif

        assert(gkp_wb<T>(0).size() == gkp_x<T>(0).size()+1);
        T area1 = gkp_wb<T>(0).back() * f_center;
        std::vector<T> fv1, fv2;
        fv1.reserve(2*gkp_x<T>(0).size()+1);
        fv2.reserve(2*gkp_x<T>(0).size()+1);
        for (uint k=0; k<gkp_x<T>(0).size(); k++) {
            const T abscissa = half_length * gkp_x<T>(0)[k];
            const T fval1 = func(center - abscissa);
            const T fval2 = func(center + abscissa);
            area1 += gkp_wb<T>(0)[k] * (fval1+fval2);
            fv1.push_back(fval1);
            fv2.push_back(fval2);
            if (retfxmap) {
                (*retfxmap)[center-abscissa] = fval1;
                (*retfxmap)[center+abscissa] = fval2;
            }
        }
#ifdef COUNTFEVAL
        nfeval+=gkp_x<T>(0).size()*2;
#endif

        integ_dbg2<<"level 0 rule: area = "<<area1<<std::endl;

        T err=0; 
        bool calcabsasc = true;
        T resabs=0., resasc=0.;
        for (int level=1; level<NGKPLEVELS; level++) {
            assert(gkp_wa<T>(level).size() == fv1.size());
            assert(gkp_wa<T>(level).size() == fv2.size());
            assert(gkp_wb<T>(level).size() == gkp_x<T>(level).size()+1);
            T area2 = gkp_wb<T>(level).back() * f_center;
            // resabs = approximation to integral of abs(f)
            if (calcabsasc) resabs = fabs(area2);
            for (uint k=0; k<fv1.size(); k++) {
                area2 += gkp_wa<T>(level)[k] * (fv1[k]+fv2[k]);
                if (calcabsasc) resabs += 
                    gkp_wa<T>(level)[k] * (fabs(fv1[k]) + fabs(fv2[k]));
            }
            for (uint k=0; k<gkp_x<T>(level).size(); k++) {
                const T abscissa = half_length * gkp_x<T>(level)[k];
                const T fval1 = func(center - abscissa);
                const T fval2 = func(center + abscissa);
                const T fval = fval1 + fval2;
                area2 += gkp_wb<T>(level)[k] * fval;
                if (calcabsasc) resabs += 
                    gkp_wb<T>(level)[k] * (fabs(fval1) + fabs(fval2));
                fv1.push_back(fval1);
                fv2.push_back(fval2);
                if (retfxmap) {
                    (*retfxmap)[center-abscissa] = fval1;
                    (*retfxmap)[center+abscissa] = fval2;
                }
            }
#ifdef COUNTFEVAL
            nfeval+=gkp_x<T>(level).size()*2;
#endif
            if (calcabsasc) {
                const T mean = area1*T(0.5);
                // resasc = approximation to the integral of abs(f-mean) 
                resasc = gkp_wb<T>(level).back() * fabs(f_center-mean);
                for (uint k=0; k<gkp_wa<T>(level).size(); k++) {
                    resasc += gkp_wa<T>(level)[k] *
                        (fabs(fv1[k]-mean) + fabs(fv2[k]-mean));
                }
                for (uint k=0; k<gkp_x<T>(level).size(); k++) {
                    resasc += gkp_wb<T>(level)[k] *
                        (fabs(fv1[k]-mean) + fabs(fv2[k]-mean));
                }
                resasc *= abs_half_length ;
                resabs *= abs_half_length;
            }
            area2 *= half_length;
            err = rescale_error (fabs(area2-area1), resabs, resasc) ;
            if (err < resasc) calcabsasc = false;

            integ_dbg2<<"at level "<<level<<" area2 = "<<area2<<
                " +- "<<err<<std::endl;

            //   test for convergence.
            if (err < epsabs || err < epsrel * fabs (area2)) {
                reg.SetArea(area2,err);
                return true;
            }
            area1 = area2;
        }

        // failed to converge 
        reg.SetArea(area1,err);

        integ_dbg2<<"failed to reach tolerance with highest-order GKP rule\n";
        return false;
    }

    template <class T>
    T FindPeriod(const std::vector<T>& v) 
    {
        uint minn = v.size()/2;
        uint n=1;
        while (n<=minn) n*=2;

        std::vector<T> dv(v.begin() + v.size()-n,v.end());
        assert(dv.size() == n);
        // Welch window:
        uint no2 = n/2;
        for(uint i=0;i<no2;i++) {
            T w = (no2-i)/no2;
            w = T(1)-w*w;
            dv[i] *= w;
            if (i>0) dv[n-i] *= w;
        }
        std::vector<std::complex<T> > fft;
        FFT(dv,fft);
        T maxpower = 0.;
        uint imax=0;
        for(uint i=1;i<fft.size();i++) {
            // start at i=1 to ignore DC power
            T power = norm(fft[i]);
            if (power > maxpower) {
                imax = i; maxpower = power;
            }
        }
        assert(imax>0);
        if (imax == 1) return 0.;

        T dimax=0., sumpower=0.;
        assert(imax >= 2);
        for(uint i=imax-2;i<=imax+2;i++) {
            if (i >= fft.size()) continue;
            T power = norm(fft[i]);
            dimax += power*i;
            sumpower += power;
        }
        assert(sumpower > 0.);
        dimax /= sumpower;
        assert(dimax > 0.);
        return T(n)/dimax;
    }

    // Slightly different if T is complex:
    template <class T>
    T FindPeriod(const std::vector<std::complex<T> >& v) 
    {
        uint minn = v.size()/2;
        uint n=1;
        while (n<=minn) n*=2;

        std::vector<std::complex<T> > fft(n);
        for(uint i=0;i<n;i++) {
            const std::complex<T>& c = v[i+v.size()-n];
            fft[i] = std::complex<T>(real(c),imag(c));
        }
        // Welch window:
        uint no2 = n/2;
        for(uint i=0;i<no2;i++) {
            T w = (no2-i)/no2;
            w = T(1)-w*w;
            fft[i] *= w;
            if(i>0) fft[n-i] *= w;
        }
        FFT(fft);
        T maxpower = 0.;
        uint imax=0;
        std::vector<T> power(n/2+2,0.);
        for(uint i=1;i<=n/2;i++) {
            // start at i=1 to ignore DC power
            power[i] = norm(fft[i])+norm(fft[n-i]);
            if (power[i] > maxpower) {
                imax = i; maxpower = power[i];
            }
        }
        if (imax == 1) return 0.;

        T dimax=0., sumpower=0.;
        for(uint i=imax-2;i<=imax+2;i++) {
            dimax += power[i]*i;
            sumpower += power[i];
        }
        dimax /= sumpower;
        return T(n)/dimax;
    }

    template <class UF>
    bool IntOsc (
        const UF& func, IntRegion<typename UF::result_type>& reg,
        const typename UF::result_type& period,
        const typename UF::result_type epsabs,
        const typename UF::result_type epsrel, 
        const typename UF::result_type extraprec=100.)
    {
        // An integration algorithm appropriate for integrals of oscillating
        // functions from a to infinity.  
        // It integrates the first 20 periods worth with GKPOSC
        // Then it successively finds the integral for one wavelength
        // at a time and extrapolates the sequence of the sum of these 
        // integrals to get the overall sum.
        // At least 10 wavelengths are calculated before trying to extrapolate
        //
        // Returns: true if succeed, false if fail
        // The area and estimated error are returned as reg.Area() and reg.Err()
        //
        typedef typename UF::result_type T;
        assert(reg.Alg() == OSCINFTWOPI || reg.Alg() == OSCINFONE);

        assert(epsabs >= 0.);
        assert(epsrel > 0.);
        if (extraprec > 2.e10) {
            std::cerr<<"Extra precision > 10^10 asked for in IntOsc.  Give up.\n";
            return false;
        }

        T x0 = reg.Left();
        T x1 = x0 + 5.*period;
        IntRegion<T> startreg(x0,x1,GKPOSC);
        if (!IntGKP(func,startreg,epsabs/extraprec,epsrel/extraprec)) {
            std::cerr<<"Integration of startreg failed in IntOsc\n";
            return false;
        }
        integ_dbg2<<"area of startreg = "<<startreg.Area()<<
            " +- "<<startreg.Err()<<std::endl;

        bool converged = false;
        const int MAXITER=64;
        SeriesLevin<T> areaseries;
        areaseries.AddTerm(startreg.Area(),startreg.Err());
        T interr = startreg.Err();
        T extraparea = startreg.Area();
        T extraperr = startreg.Err();
        for(int iter=1; iter<=MAXITER; iter++) {
            x0 = x1;
            x1 += period;
            IntRegion<T> oneperiod(x0,x1,GKPOSC);
            if (!IntGKP(func,oneperiod,epsabs/extraprec,epsrel/extraprec)) {
                std::cerr<<"Integration of oneperiod failed in IntOsc\n";
                return false;
            }
            integ_dbg2<<"area of period "<<iter<<" = "<<oneperiod.Area()<<
                " +- "<<oneperiod.Err()<<std::endl;
            areaseries.AddTerm(oneperiod.Area(),oneperiod.Err());
            interr += oneperiod.Err();
            areaseries.Extrapolate(&extraparea,&extraperr,0);
            extraperr *= 100.;  // It can be optimistic at times, up to 2 digits
            if (iter >= 5) {
                T toler = std::max(epsabs,epsrel*fabs(extraparea));
                integ_dbg2<<"area after "<<iter<<" oscillations = "<<extraparea;
                integ_dbg2<<" +- "<<extraperr<<" (tol = "<<toler<<")\n";
                if (extraperr < toler) {
                    converged = true;
                    break;
                }
                if (interr > toler) {
                    integ_dbg2<<"interr = "<<interr<<", epsabs = "<<epsabs<<
                        ", epsrel*area = "<<epsrel*fabs(extraparea)<<"\n"
                        "Integration errors already too big to ever converge\n"
                        "Redo with extraprec = "<<extraprec*10<<std::endl;
                    return IntOsc(func,reg,period,epsabs,epsrel,extraprec*10);
                }
                if (areaseries.PastBest()) {
                    integ_dbg2<<
                        "Series extrapolation is not going to improve.\n"
                        "Redo with extraprec = "<<extraprec*10<<std::endl;
                    return IntOsc(func,reg,period,epsabs,epsrel,extraprec*10);
                }
            } else {
                integ_dbg2<<"area after "<<iter<<" oscillations = "<<extraparea;
                integ_dbg2<<" +- "<<extraperr<<std::endl;
            }
        }
        reg.SetArea(extraparea,extraperr);
        if (!converged) {
            integ_dbg2<<
                "Series extrapolation not good enough after MAXITER "
                "iterations\n";
            // Probably the period is not the best period to use
            // Find a better (longer) period for another pass:
            T newperiod = FindPeriod(areaseries.GetTerms());
            integ_dbg2<<"newperiod = "<<newperiod<<std::endl;
            if (newperiod == 0.0) {
                integ_dbg2<<
                    "There isn't a better period for this.  "
                    "This one is just a kicker.\n";
                newperiod = 100.;
            }
            integ_dbg2<<"Redo with period = "<<newperiod<<
                "*period (= "<<newperiod*period<<")\n";
            return IntOsc(func,reg,newperiod*period,epsabs,epsrel,extraprec);
        }
        return true;
    }

    template <class UF>
    bool IntGKP (const UF& func, IntRegion<typename UF::result_type>& reg,
                 const typename UF::result_type epsabs,
                 const typename UF::result_type epsrel,
                 std::map<typename UF::result_type,
                 typename UF::result_type>* retfxmap=0)
    {
        // An adaptive integration algorithm which computes the integral of f
        // over the region reg.
        // First the non-adaptive GKP algorithm is tried.
        // If that is not accurate enough (according to the absolute and
        // relative accuracies, epsabs and epsrel),
        // the region is split in half, and each new region is integrated.
        // The routine continues by successively splitting the subregion
        // which gave the largest absolute error until the integral converges.
        //
        // Returns true is succeed, false if fail
        // The area and estimated error are returned as reg.Area() and reg.Err()
        // If desired, *retx and *retf return std::vectors of x,f(x) 
        // respectively
        // They only include the evaluations in the non-adaptive pass, so they
        // do not give an accurate estimate of the number of function 
        // evaluations.
        typedef typename UF::result_type T;
        bool oscsplit = (reg.Alg() == GKPOSC);
        integ_dbg2<<"Start IntGKP with oscsplit = "<<oscsplit<<std::endl;
        assert(oscsplit || reg.Alg() == GKP);

        assert(epsabs >= 0.);
        assert(epsrel > 0.);

        // perform the first integration 
        bool done;
        std::map<T,T> fxmap;
        if (oscsplit) {
            done = IntGKPNA(func, reg, epsabs, epsrel, &fxmap);
            if (!done) reg.FindOscSplits(fxmap);
        } else {
            done = IntGKPNA(func, reg, epsabs, epsrel, &fxmap);
        }
        if (retfxmap) *retfxmap = fxmap;
        if (done) return true;

        integ_dbg2<<"In adaptive GKP, failed first pass... subdividing\n";
        integ_dbg2<<"Intial range = "<<reg.Left()<<".."<<reg.Right()<<std::endl;

        int roundoff_type1 = 0, error_type = 0;
        T roundoff_type2 = 0.;
        uint iteration = 1;

        std::priority_queue<IntRegion<T>,std::vector<IntRegion<T> > > allregions;
        allregions.push(reg);
        T finalarea = reg.Area();
        T finalerr = reg.Err();
        T tolerance= std::max(epsabs, epsrel * fabs(finalarea));
        assert(finalerr > tolerance);

        while(!error_type && finalerr > tolerance) {
            // Bisect the subinterval with the largest error estimate 
            integ_dbg2<<"Current answer = "<<finalarea<<
                " +- "<<finalerr<<"  (tol = "<<tolerance<<")\n";
            IntRegion<T> parent = allregions.top(); 
            allregions.pop();
            integ_dbg2<<"Subdividing largest error region "<<parent.Left()<<
                ".."<<parent.Right()<<std::endl;
            integ_dbg2<<"parent area = "<<parent.Area()<<
                " +- "<<parent.Err()<<std::endl;
            std::vector<IntRegion<T> > children;
            parent.SubDivide(&children);
            // For "GKP", there are only two, but for GKPOSC, there is one 
            // for each oscillation in region

            // Try to do at least 3x better with the children
            T factor = 3*children.size()*finalerr/tolerance;
            T newepsabs = fabs(parent.Err()/factor);
            T newepsrel = newepsabs/fabs(parent.Area());
            integ_dbg2<<"New epsabs,rel = "<<newepsabs<<','<<newepsrel<<
                "  ("<<children.size()<<" children)\n";

            T newarea = T(0.0);
            T newerror = 0.0;
            for(uint i=0;i<children.size();i++) {
                IntRegion<T>& child = children[i];
                integ_dbg2<<"Integrating child "<<child.Left()<<
                    ".."<<child.Right()<<std::endl;
                bool converged;
                if (oscsplit) {
                    converged = IntGKPNA(
                        func, child, newepsabs, newepsrel, &fxmap);
                    if (!converged) child.FindOscSplits(fxmap);
                } else {
                    converged = IntGKPNA(func, child, newepsabs, newepsrel);
                }
                integ_dbg2<<"child ("<<i+1<<'/'<<children.size()<<") ";
                if (converged) {
                    integ_dbg2<<" converged.";
                } else {
                    integ_dbg2<<" failed.";
                }
                integ_dbg2<<"  Area = "<<child.Area()<<
                    " +- "<<child.Err()<<std::endl;

                newarea += child.Area();
                newerror += child.Err();
            }
            integ_dbg2<<"Compare: newerr = "<<newerror<<
                " to parent err = "<<parent.Err()<<std::endl;

            finalerr += (newerror - parent.Err());
            finalarea += newarea - parent.Area();

            T delta = parent.Area() - newarea;
            if (newerror <= parent.Err() && fabs (delta) <=  parent.Err() 
                && newerror >= 0.99 * parent.Err()) {
                integ_dbg2<<
                    "roundoff type 1: delta/newarea = "<<
                    fabs(delta)/fabs(newarea)<<
                    ", newerror/error = "<<newerror/parent.Err()<<std::endl;
                roundoff_type1++;
            }
            if (iteration >= 10 && newerror > parent.Err() 
                && fabs(delta) <= newerror-parent.Err()) {
                integ_dbg2<<"roundoff type 2: newerror/error = "<<
                    newerror/parent.Err()<<std::endl;
                roundoff_type2+=std::min(newerror/parent.Err()-1.,T(1.));
            }

            tolerance = std::max(epsabs, epsrel * fabs(finalarea));
            if (finalerr > tolerance) {
                if (roundoff_type1 >= 200) {
                    error_type = 1;	// round off error 
                    integ_dbg2<<"GKP: Round off error 1\n";
                }
                if (roundoff_type2 >= 200.) {
                    error_type = 2;	// round off error 
                    integ_dbg2<<"GKP: Round off error 2\n";
                }
                const double parent_size = parent.Right()-parent.Left();
                const double reg_size = reg.Right()-parent.Left();
                if (fabs(parent_size / reg_size) < Epsilon<double>()) {
                    error_type = 3; // found singularity
                    integ_dbg2<<"GKP: Probable singularity\n";
                }
            }
            for(uint i=0;i<children.size();i++) allregions.push(children[i]);
            iteration++;
        } 

        // Recalculate finalarea in case there are any slight rounding errors
        finalarea=0.; finalerr=0.;
        while (!allregions.empty()) {
            const IntRegion<T>& r=allregions.top();
            finalarea += r.Area();
            finalerr += r.Err();
            allregions.pop();
        }
        reg.SetArea(finalarea,finalerr);

        if (error_type == 1) {
            std::cerr << 
                "Type 1 roundoff's = "<<roundoff_type1<<
                ", Type 2 = "<<roundoff_type2<<"\n"
                "Roundoff error 1 prevents tolerance from being achieved in "
                "IntGKP\n";
            return false;
        } else if (error_type == 2) {
            std::cerr << 
                "Type 1 roundoff's = "<<roundoff_type1<<
                ", Type 2 = "<<roundoff_type2<<"\n"
                "Roundoff error 2 prevents tolerance from being achieved in "
                "IntGKP\n";
            return false;
        } else if (error_type == 3) {
            std::cerr << 
                "Bad integrand behavior found in the integration interval in "
                "IntGKP\n";
            return false; 
        } else {
            return true;
        }
    }


    const double DEFABSERR = 1.e-15;
    const double DEFRELERR = 1.e-6;

    template <class UF>
    struct AuxFunc1 :   // f(1/x-1) for int(a..infinity)
        public std::unary_function<typename UF::argument_type,
        typename UF::result_type> 
    {
    public:
        AuxFunc1(const UF& _f) : f(_f) {}
        typename UF::result_type 
            operator()(typename UF::argument_type x) const 
            { return f(1./x-1.)/(x*x); }
    private:
        const UF& f;
    };

    template <class UF>
    struct AuxFunc1<UF> Aux1(UF uf) 
    { return AuxFunc1<UF>(uf); }

    template <class UF>
    struct AuxFunc2 :  // f(1/x+1) for int(-infinity..b)
        public std::unary_function<typename UF::argument_type,
        typename UF::result_type> 
    {
    public:
        AuxFunc2(const UF& _f) : f(_f) {}
        typename UF::result_type 
            operator()(typename UF::argument_type x) const 
            { return f(1./x+1.)/(x*x); }
    private:
        const UF& f;
    };

    template <class UF>
    struct AuxFunc2<UF> Aux2(UF uf) 
    { return AuxFunc2<UF>(uf); }

    template <class UF>
    struct AuxFunc3 :  // f(-x) for int(-infinity..b) for OSCINF
        public std::unary_function<typename UF::argument_type,
        typename UF::result_type> 
    {
    public:
        AuxFunc3(const UF& _f) : f(_f) {}
        typename UF::result_type 
            operator()(typename UF::argument_type x) const 
            { return f(-x); }
    private:
        const UF& f;
    };

    template <class UF>
    struct AuxFunc3<UF> Aux3(UF uf) 
    { return AuxFunc3<UF>(uf); }

    template <class UF>
    typename UF::result_type int1d(
        const UF& func, 
        IntRegion<typename UF::result_type>& reg,
        const typename UF::result_type& abserr=DEFABSERR,
        const typename UF::result_type& relerr=DEFRELERR)
    {
        typedef typename UF::result_type T;

        integ_dbg2<<"start int1d: "<<reg.Left()<<".."<<reg.Right()<<std::endl;

        const T MOCKINF = 1.e10;
        const double PI = 3.141592653589793;
        const double TWOPI = 2.*PI;

        if ((reg.Left() <= -MOCKINF && reg.Right() > 0) ||
            (reg.Right() >= MOCKINF && reg.Left() < 0)) { 
            reg.AddSplit(0);
        }

        if (reg.NSplit() > 0) {
            std::vector<IntRegion<T> > children;
            reg.SubDivide(&children);
            integ_dbg2<<"Subdivided into "<<children.size()<<" children\n";
            T answer=T();
            T err; 
            for(uint i=0;i<children.size();i++) {
                IntRegion<T>& child = children[i];
                integ_dbg2<<"i = "<<i<<": bounds = "<<child.Left()<<
                    ','<<child.Right()<<std::endl;
                answer += int1d(func,child,abserr,relerr);
                err += child.Err();
                integ_dbg2<<"subint = "<<child.Area()<<
                    " +- "<<child.Err()<<std::endl;
            }
            reg.SetArea(answer,err);
            return answer;
        } else {
            bool fail=false;
            T period=0.;
            switch(reg.Alg()) {
                // Any additional OSCINF periods can go here and fall through
              case OSCINFTWOPI : {
                  if(reg.Alg() == OSCINFTWOPI) period = TWOPI;
              } // fall through...
              case OSCINFONE : {
                  if(reg.Alg() == OSCINFONE) period = 1.;
                  integ_dbg2<<"oscinf, period = "<<period<<std::endl;
                  if (reg.Left() <= -MOCKINF) {
                      integ_dbg2<<"left = -infinity, right = "<<
                          reg.Right()<<std::endl;
                      assert(reg.Right() <= 0.);
                      IntRegion<T> modreg(
                          -reg.Right(),-reg.Left(),reg.Alg(),reg.dbgout);
                      if (!IntOsc(Aux3(func),modreg,period,abserr,relerr)) { 
                          fail = true; std::cerr<<"IntOsc failed\n"; 
                      }
                      reg.SetArea(modreg.Area(),modreg.Err());
                  } else if (reg.Right() >= MOCKINF) {
                      integ_dbg2<<"left = "<<reg.Left()<<", right = infinity\n";
                      assert(reg.Left() >= 0.);
                      if(!IntOsc(func,reg,period,abserr,relerr))
                      { fail = true; std::cerr<<"IntOsc failed\n"; }
                  } else {
                      integ_dbg2<<"left = "<<reg.Left()<<", right = "<<
                          reg.Right()<<std::endl;
                      IntRegion<T> modreg(reg.Left(),reg.Right(),GKPOSC,
                                          reg.dbgout);
                      if (!IntGKP(func,modreg,abserr,relerr)) { 
                          fail = true; std::cerr<<"IntGKP failed\n"; 
                      }
                      reg.SetArea(modreg.Area(),modreg.Err());
                  }
              } break;
              case GKP : case GKPOSC : {
                  if (reg.Left() <= -MOCKINF) {
                      integ_dbg2<<"left = -infinity, right = "<<
                          reg.Right()<<std::endl;
                      assert(reg.Right() <= 0.);
                      IntRegion<T> modreg(1./(reg.Right()-1.),0.,reg.Alg(),
                                          reg.dbgout);
                      if (!IntGKP(Aux2<UF>(func),modreg,abserr,relerr)) {
                          fail = true; std::cerr<<"IntGKP failed\n"; 
                      }
                      reg.SetArea(modreg.Area(),modreg.Err());
                  } else if (reg.Right() >= MOCKINF) {
                      integ_dbg2<<
                          "left = "<<reg.Left()<<", right = infinity\n";
                      assert(reg.Left() >= 0.);
                      IntRegion<T> modreg(
                          0.,1./(reg.Left()+1.),reg.Alg(),reg.dbgout);
                      if(!IntGKP(Aux1<UF>(func),modreg,abserr,relerr)) {
                          fail = true; std::cerr<<"IntGKP failed\n"; 
                      }
                      reg.SetArea(modreg.Area(),modreg.Err());
                  } else {
                      integ_dbg2<<"left = "<<reg.Left()<<
                          ", right = "<<reg.Right()<<std::endl;
                      if (!IntGKP(func,reg,abserr,relerr)) {
                          fail = true; std::cerr<<"IntGKP failed\n"; 
                      }
                  }
              } break;
              default : {
                  myerror("Invalid algorithm");
              } break;
            }
            if (fail) myerror("Integration failed"); 
            integ_dbg2<<"done int1d  answer = "<<reg.Area()<<
                " +- "<<reg.Err()<<std::endl;
            return reg.Area();
        }
    }

    template <class BF, class YREG>
    class Int2DAuxType : 
        public std::unary_function<typename BF::first_argument_type,
        typename BF::result_type> 
    {
    public:
        Int2DAuxType(const BF& _func, const YREG& _yreg,
                     const typename BF::result_type& _abserr,
                     const typename BF::result_type& _relerr) :
            func(_func),yreg(_yreg),abserr(_abserr),relerr(_relerr) 
        {}

        typename BF::result_type operator()(
            typename BF::first_argument_type x) const 
        {
            typename YREG::result_type tempreg = yreg(x);
            typename BF::result_type result = 
                int1d(bind21(func,x),tempreg,abserr,relerr);
            integ_dbg3<<"Evaluated int2dAux at x = "<<x<<": f = "<<result<<
                " +- "<<tempreg.Err()<<std::endl;
            return result;
        } 

    private:
        const BF& func;
        const YREG& yreg;
        typename BF::result_type abserr,relerr;
    };

    template <class BF, class YREG>
    typename BF::result_type int2d(
        const BF& func,
        IntRegion<typename BF::result_type>& reg,
        const YREG& yreg, const typename BF::result_type& abserr=DEFABSERR,
        const typename BF::result_type& relerr=DEFRELERR)
    {
        integ_dbg2<<"Starting int2d: range = "<<reg.Left()<<".."<<
            reg.Right()<<std::endl;
        Int2DAuxType<BF,YREG> faux(func,yreg,abserr*1.e-3,relerr*1.e-3);
        typename BF::result_type answer = int1d(faux,reg,abserr,relerr);
        integ_dbg2<<"done int2d  answer = "<<answer<<
            " +- "<<reg.Err()<<std::endl;
        return answer;
    }

    template <class TF, class YREG, class ZREG>
    class Int3DAuxType : 
        public std::unary_function<typename TF::firstof3_argument_type,
        typename TF::result_type> 
    {
    public:
        Int3DAuxType(
            const TF& _func, const YREG& _yreg, const ZREG& _zreg, 
            const typename TF::result_type& _abserr,
            const typename TF::result_type& _relerr) :
            func(_func),yreg(_yreg),zreg(_zreg),
            abserr(_abserr),relerr(_relerr) 
        {}

        typename TF::result_type operator()(
            typename TF::firstof3_argument_type x) const 
        {
            typename YREG::result_type tempreg = yreg(x);
            typename TF::result_type result = 
                int2d(bind31(func,x),tempreg,bind21(zreg,x),abserr,relerr);
            integ_dbg3<<"Evaluated int3dAux at x = "<<x<<
                ": f = "<<result<<" +- "<<tempreg.Err()<<std::endl;
            return result;
        }

    private:
        const TF& func;
        const YREG& yreg;
        const ZREG& zreg;
        typename TF::result_type abserr,relerr;
    };

    template <class TF, class YREG, class ZREG>
    typename TF::result_type int3d(
        const TF& func,
        IntRegion<typename TF::result_type>& reg,
        const YREG& yreg, const ZREG& zreg,
        const typename TF::result_type& abserr=DEFABSERR,
        const typename TF::result_type& relerr=DEFRELERR)
    {
        integ_dbg2<<"Starting int3d: range = "<<reg.Left()<<".."<<
            reg.Right()<<std::endl;
        Int3DAuxType<TF,YREG,ZREG> faux(
            func,yreg,zreg,abserr*1.e-3,relerr*1.e-3);
        typename TF::result_type answer = int1d(faux,reg,abserr,relerr);
        integ_dbg2<<"done int3d  answer = "<<answer<<
            " +- "<<reg.Err()<<std::endl;
        return answer;
    }

    // Helpers for constant regions for int2d, int3d:

    template <class T>
    struct ConstantReg1 : std::unary_function<T, IntRegion<T> >
    {
        ConstantReg1(T a,T b,INTALG alg=DEFALG) : ir(a,b,alg) {}
        IntRegion<T> operator()(T x) const {return ir;}
        IntRegion<T> ir;
    };

    template <class T>
    struct ConstantReg2 : std::binary_function<T, T, IntRegion<T> >
    {
        ConstantReg2(T a,T b,INTALG alg=DEFALG) : ir(a,b,alg) {}
        IntRegion<T> operator()(T x, T y) const {return ir;}
        IntRegion<T> ir;
    };

}

#endif
